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ABSTRACT 

Analysis  of  the  Effect  of  Matrix  Degradation  on  Fatigue 
Behavior  of  a  Graphite/ Epoxy  Laminate.  (May  1982) 

Robert  Thomas  Arenburg,  B.C.E.T.,  Southern  Technical  Institute 
Chairman  of  Advisory  Committee:  Dr.  R.  A.  Schapery 

An  analytical  study  of  the  effect  of  fatigue  on  a  [±45/902]$ 
graphite/epoxy  laminate  is  conducted  considering  two  levels  of  damage. 
The  first  one  is  microscopic  damage,  which  is  viewed  as  a  matrix  soft¬ 
ening  mechanism  that  is  accounted  for  by  the  use  of  micromechanics. 

The  second  level  of  damage  is  delamination.  An  evaluation  of  lamina¬ 
tion  theory  versus  a  finite  element  elasticity  solution  which  accounts 
for  free  edge  stresses  is  made  first.  Lamination  theory  is  found  to 
be  inadequate  due  to  its  inability  to  account  for  large  amounts  of 
shear  deformation  as  the  matrix  degrades.  Experimental  data  are  used 
to  evaluate  the  matrix  degradation  as  a  function  of  fatigue  cycles, 
and  then  a  parametric  study  of  the  strain  energy  release  rate  as  a 
function  of  delamination  size  at  various  degrees  of  matrix  degradation 
is  conducted.  For  a  constant  maximum  stress  fatigue  test,  the 
matrix  degradation  leads  to  an  increase  in  strain  energy  which  is  pre¬ 
dicted  to  eventually  cause  rapid  delamination.  Reasonably  good  agree¬ 
ment  between  theory  and  experimental  data  on  the  [±45/902^  laminate 
is  demonstrated. 
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I.  INTRODUCTION 

Fiber  reinforced  composite  materials  are  known  for  their  good 
fatigue  characteristics.  Helicopter  rotor  blade  technology  has  recog 
nized  this  fact,  and  today  most  current  rotor  blade  designs  have  adopt 
ed  the  use  of  composite  materials.  The  study  conducted  considers  the 
influence  of  two  fatigue  mechanisms:  matrix  degradation  and  free  edge 
delaminations,  and  their  influence  on  deformation  and  failure  behavior. 

Consider  a  matrix  dominated,  graphite/epoxy  laminate,  which  is 
fatigued  under  constant  maximum  stress  conditions.  The  mean  strain 
response  as  a  function  of  fatigue  cycles  has  been  characterized  as 
predominately  linear  if  plotted  on  log-log  scales.  Near  failure  the 
mean  strain  departs  from  its  previous  linear  path  and  increases  very 
rapidly  to  failure.  It  is  speculated  that  the  laminate  response  to 
the  fatigue  loading  can  be  explained  as  follows. 

The  initial  linear  response  is  a  result  of  an  accumulation  of 
matrix  "damage."  The  reason  for  the  damage  to  be  associated  with  the 
matrix  can  be  explained  by  the  considerably  higher  stiffness  and 
strength  of  the  graphite  fibers  compared  to  the  matrix.  When  consid¬ 
ering  a  matrix  dominated  laminate,  i.e.,  a  laminate  without  fibers 
oriented  in  the  direction  of  the  load  introduction,  the  fibers  can 
be  viewed  as  rigid  inclusions  and  thus  act  as  stress  concentrators 
or  stress  risers.  Due  to  the  stress  riser  effect  of  the  fibers  and 
their  higher  strength,  the  matrix  is  the  first  phase  damaged  by  the 
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fatigue  loading.  An  important  form  of  damage  induced  by  the  rigid 
inclusion  effect  would  be  a  fiber-matrix  disbond. 

Polymers  with  a  highly  nonuniform  crosslink  density  could  pro¬ 
vide  "soft  spots,"  weak  regions  causing  strain  concentrations.  The 
soft  spots  could  provide  initiation  sites  for  microcracking.  Other 
zones  of  damage  could  be  composed  of  microscopic  voids  and  cracks. 

In  general,  the  term  "damage"  will  be  used  in  reference  to  micro¬ 
scopic  events  which  result  in  reduced  global  laminate  stiffness. 

The  next  phase  might  be  the  interaction  of  the  matrix  damage, 
characterized  by  microscopic  flaws,  with  macroscopic  delaminations. 
Once  the  delamination  has  penetrated  most  of  the  specimen  width,  fail¬ 
ure  should  occur  very  rapidly  due  to  the  degraded  matrix  strength. 

The  prime  concern  of  this  thesis  is  to  account  for  the  observed 
fatigue  behavior  through  an  analytical  study  of  the  effects  of  matrix 
degradation  and  free  edge  delaminations. 

Present  Status 

Frequently  the  fracture  damage  in  fibrous  composite  laminates 
involving  matrix  cracks  can  be  characterized  by  two  damage  modes. 

The  first  mode  consists  of  multiple  transverse  cracks  which  occur  in 
plies  with  fibers  oriented  perpendicular  or  nearly  perpendicular  to 
the  applied  load.  The  second  mode  consists  of  inter-ply  delamina¬ 
tions  which  usually  start  near  the  free  edge  of  the  laminate. 

Experimental  studies  have  shown  there  exists  a  threshold  stress 
state  below  which  neither  transverse  cracks  nor  edge  delaminations 
occur.  Bader  et  al.  [1]  showed  that  the  strain  required  for  initia- 
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tion  of  transverse  cracks  in  a  [0n/90m]s  graphite/epoxy  laminate  in¬ 
creased  as  the  thickness  of  the  90°  plies  decreased.  The  outside  0° 
plies  of  the  laminate  performed  as  an  arresting  barrier  for  the  trans¬ 
verse  cracks.  Also  noted  was  a  state  of  crack  saturation  which  was 
characterized  by  a  typical  crack  spacing.  These  results  confirm  the 
idea  of  ply  interaction  in  the  transverse  cracking  process.  The  re¬ 
sults  of  Bader  et  al.  confirm  the  earlier  trends  noted  in  a  microme¬ 
chanical  analysis  by  Aveston  and  Kelly  [2]. 

Free  edge  delamination  has  been  attributed  to  locally  high  edge 
stresses,  in  which  the  magnitude  and  direction  of  the  stresses  are 
dependent  on  the  ply  stacking  sequence.  Rodine  and  Eisenmann  [3] 
showed  for  a  C±45n/0n/90n]$  graphite/epoxy  laminate  that  the  onset  of 
edge  delamination  occurred  at  a  lower  stress  level  for  thicker  lami¬ 
nates.  They  used  a  statistical  argument  that  a  thicker  laminate  had 
potentially  more  flaws  to  initiate  the  delamination. 

We  believe  the  observed  thickness  effect  on  the  threshold  stress 
for  edge  delamination  and  transverse  cracking  may  be  explained  from 
an  energy  argument.  The  strain  energy  contained  within  a  laminate 
plays  an  important  role  in  the  initiation  and  behavior  of  cracks  and 
from  this  the  ply  thickness  becomes  an  important  parameter.  Since 
the  strain  energy  is  dependent  upon  the  ply  interaction,  and  this 
can  be  altered  by  the  kinematics  of  a  crack,  the  study  of  the  energy 
state  of  a  laminate  should  provide  a  criterion  for  crack  growth. 

Rybicki  et  al.  [4]  used  linear  fracture  mechanics  in  character¬ 
izing  free  edge  delamination.  The  critical  strain  energy  release  rate 
was  used  as  a  criterion  for  unstable  crack  growth.  In  their  analysis. 
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free  edge  flaws  were  assumed  to  initiate  the  delamination.  The  finite 
element  method  in  conjunction  with  a  virtual  crack  closure  method  [5] 
was  used  to  calculate  the  energy  release  rate.  The  energy  release 
rate  was  found  to  be  constant  during  stable  crack  growth,  and  was 
based  upon  extensional  strain  and  delamination  size  provided  by  uni¬ 
axial  tensile  tests.  From  these  results  it  was  verified  that  linear 
elastic  fracture  mechanics  could  be  applied  to  fiber  reinforced  com¬ 
posites  and  the  strain  energy  release  rate  could  be  a  suitable  param¬ 
eter. 

Wang  and  Crossman  et.  al .  [6,  7]  continued  the  energy  approach, 
recognizing  its  ability  to  account  for  the  structural  interaction  of 
plies.  They  considered  both  transverse  cracking  and  edge  delamination 
in  uniaxial  tensile  specimens.  Energy  release  rate,  G,  as  a  function 
of  crack  length  due  to  a  unit  load  was  introduced.  This  function  was 
used  to  describe  the  energy  available  for  crack  growth.  The  finite 
element  method  was  used  to  evaluate  the  function;  from  this  the  thresh¬ 
old  stresses  for  edge  delamination  and  transverse  cracking  were  calcu¬ 
lated  for  a  constant  critical  G,  which  yielded  the  ply  thickness  depen¬ 
dence  of  the  threshold  stress  for  transverse  cracking  and  edge  delami¬ 
nation.  Correlation  of  their  method  with  experimental  data  of  a 
[±25/90nJs  (n  =  1,2,3)  graph ite/epoxy  laminate  showed  good  agreement. 

O'Brien  [8]  developed  a  simple  energy  approach  to  analyze  delami¬ 
nation  onset  and  growth.  His  analysis  used  classical  lamination  the¬ 
ory  to  evaluate  the  strain  energy  release  rate  based  upon  a  simple 
rule-of-mixtures  expression  for  the  laminate  stiffness.  In  studying 
a  [±30/±30/90/90sJs  laminate,  O'Brien  did  not  try  to  model  both  trans- 
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verse  cracking  in  the  90°  plies  and  interlaminar  delaminations  explic¬ 
itly.  Instead,  he  considered  an  idealized  delamination,  as  shown  in 
Figure  1.  The  modeled  delamination.  Figure  1(b),  consists  of  two 
planar  delaminations  at  the  -30/90  interfaces.  The  reason  for  using 
the  two  delaminations  in  the  model  is  to  account  for  two  mechanisms. 

The  first  mechanism  is  the  delamination  shifting  from  one  -30/90 
interface  to  the  other  through  existing  transverse  cracks  in  the  90° 
plies.  This  shifting  is  believed  to  reduce  the  bending-extension 
coupling  effect  which  would  occur  if  the  delamination  was  planar  at 
only  one  of  the  -30/90  ply  interfaces.  Hence,  the  two  symmetric 
delaminations  were  used  to  reduce  the  bending-extension  coupling 
effect.  The  second  mechanism  is  the  effect  of  the  transverse  cracks 
in  the  90°  plies.  The  transverse  cracks  reduce  the  axial  stiffness 
of  the  laminate;  this  effect  was  small  and  was  considered  "lumped" 
into  the  use  of  the  two  symmetric  delaminations. 

Micromechanics  has  been  used  to  model  the  limiting  strength  of 
composites  by  reducing  the  matrix  modulus  while  retaining  the  original 
modulus  of  the  reinforcing  fibers  [9,  10].  In  Reference  [9],  the 
static  failure  stress  of  a  short  glass-fiber  reinforced  polyester 
composite  was  modeled  by  reducing  the  matrix  modulus  to  zero.  This 
produced  results  in  close  agreement  with  experimental  data.  In 
Reference  [10],  the  large  amounts  of  strain  which  occurred  during 
fatigue  testing  was  partially  accounted  for  by  reducing  the  matrix 
modulus  to  zero.  The  micromechanics  model  accounted  for  the  mean 
strain  response  very  well  until  the  strain  increased  rapidly  close  to 
failure.  The  departure  of  the  mean  strain  from  the  micromechanics 


(b)  Modeled  Delamination. 


Figure  1  .  Idealized  model  of  observed  damage. 
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model  was  speculated  to  be  due  to  other  failure  mechanisms  such  as 
delamination. 

Method  of  Analysis 

In  this  study,  an  in-depth  examination  of  the  damage  state  of 
fatigued  composite  laminates  will  be  performed.  Since  in  Reference 
[10]  it  was  shown  that  not  all  of  the  strain  observed  could  be  account¬ 
ed  for  by  reducing  the  matrix  modulus  in  a  micromechanics  model,  it  is 
felt  that  a  more  detailed  study  of  the  same  [±45/902]2  graphite/epoxy 
laminate  is  required. 

The  first  phase  of  this  study  involves  an  elasticity  solution 
for  the  composite.  The  finite  element  method  will  be  employed  for 
this  due  to  the  complexity  of  the  problem.  A  state  of  "uniform  axial 
extension"  as  described  by  Reference  [11]  will  be  assumed.  This  is 
the  same  as  the  "generalized  plane  strain  state"  used  in  Reference 
[12]. 

The  micromechanics  equations  for  fibrous  composites  presented  in 
[13]  will  be  used  to  obtain  the  effective  lamina  moduli  as  functions 
of  the  constituent  properties  and  the  fiber  volume  fraction.  The 
matrix  modulus  will  be  reduced  uniformly  throughout  the  laminate  in 
order  to  account  for  fatigue  induced  matrix  damage.  It  has  been 
suggested  [14]  that  the  matrix  modulus  reduction  could  be  a  function 
of  the  local  average  octahedral  shear  stress  in  the  matrix.  A  future 
study  of  this  should  be  made  to  account  for  stress  gradients. 

The  second  phase  of  the  study  will  involve  a  fracture  mechanics 
analysis  of  the  laminate  with  reduced  matrix  modulus,as  it  is  doubtful 
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that  the  modulus  reduction  alone  can  account  for  the  entire  fatigue 
behavior  of  the  composite.  Also,  the  analysis  method  of  O'Brien  [8] 
will  be  used  as  a  comparison  with  the  finite  element  based  analysis. 

The  search  for  a  locally  available  finite  element  program  to  per¬ 
form  the  required  stress  and  fracture  analyses  was  quickly  abandoned 
since  the  available  programs  lacked  anisotropic  material  behavior  and 
fracture  analysis.  A  simple  two-dimensional,  plane  stress/plane  strain 
program  with  an  eight  node  isoparametric  element  listed  in  Reference 
[15]  was  chosen  for  development.  For  crack  analysis,  the  virtual 
crack  closure  method  as  described  in  Reference  [5]  is  used.  This 
method  is  particularly  attractive,  since  a  special  crack  tip  element 
is  not  required. 

The  experimental  data  used  for  correlation  is  from  the  extensive 
characterization  of  AS/3501-6  graphite/epoxy  performed  by  Dr.  T.  Ho 
of  Vought  Corporation,  Dallas,  Texas  [10,16,17]. 
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II.  MICROMECHANICS 

Micromechanics  provides  the  means  by  which  the  homogeneous  consti¬ 
tuent  properties  of  a  composite  are  related  to  the  macroscopically 

I 

observed  properties.  This  section  will  provide  the  means  by  which  the 
matrix  modulus  can  be  related  to  the  laminate  deformation  behavior. 

The  author  makes  no  attempt  to  cover  the  vast  amount  of  material  in 
the  field  of  micromechanics  other  than  to  document  the  pertinent 
items  used  in  this  investigation.  Readers  are  referred  to  References 
[18,19,  and  20]  for  additional  information  on  the  subject  of  micro¬ 
mechanics.  The  article  by  Chamis  and  Sendeckyj  [18]  is  an  excellent 
overview  of  many  analytic  micromechanics  theories.  Hashin  presents  an 
extensive  treatment  of  the  subject  in  Reference  [19].  More  recent 
work  in  micromechanics  has  begun  to  rely  upon  numerical  methods  to 
model  inelastic  behavior;  a  brief  history  of  these  are  presented  in 
Reference  [20]. 

The  Constitutive  Relationship 

The  constitutive  relationship  for  a  linearly  elastic  body  includ¬ 
ing  thermal  effects  can  be  expressed  in  contracted  notation  as: 

ai  =  Cij  ^ej  "  ajAT^  0»J  =  ]»6)  (1) 

or  inverted  as 

ei  =  s-jj  crj  +  a. AT  ( i , j  =  1,6)  (2) 

In  the  above  expressions  the  repeated  j  index  is  summed  and  i  is  a 
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free  index.  Both  indices  range  from  1  to  6.  The  quantities  o.. ,  ^ 
and  are  the  stress,  strain  and  coefficient  of  thermal  expansion 
matrices,  respectively.  These  are  defined  as: 


{a}T  = 

(axx’ 

V 

Q 

N 

N 

*0 

ayz’ 

azx5 

0xy} 

(3) 

{e)T  = 

(£xx’ 

eyy’ 

£zz’ 

eyz’ 

£zx> 

exy^ 

(4) 

{a}T  = 

(axx’ 

ayy’ 

azz* 

o  , 

0  , 

0  ) 

(5) 

The  term  AT  is  the  change  in  temperature  from  a  stress  free  reference 
temperature.  The  remaining  terms,  C..  and  S..,  are  the  modulus  matrix 

*  J  1  J 

and  compliance  matrix  components,  respectively.  These  will  be  de¬ 
fined  next. 

Consider  a  unidirectional  fiber  reinforced  laminate  as  shown  in 
Figure  2.  Notice,  for  fibers  which  are  parallel  to  axis  1  and  ran¬ 
domly  positioned,  axis  1  is  an  axis  of  isotropy.  Hence,  the  laminate 
is  transversely  isotropic  and  requires  only  five  independent  constants 
in  the  stress-strain  relationship  [21]. 

The  modulus  matrix  in  the  principal  material  coordinates  shown 
in  Figure  2  is: 


where  C44  =  1/2(C22  -  C23) 


*e  2  .  Principal  material  coordinates,  (1,2,3)  used 
to  define  the  engineering  constants. 
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Alternatively. 

,  we  can 

express  the  compl iance  matrix 

in  terms 

of  the 

engineering 

constants: 

l/E, 

"v21/E2 

_v21/E2 

0 

0 

0 

‘VEi 

i/e2 

_v32/E2 

0 

0 

0 

s..  = 

LU 

CNJ 

1 

"v23/E2 

i/e2 

0 

0 

0 

(7) 

1  J 

0 

0 

0 

1/G23 

0 

0 

0 

0 

0 

0 

1/G12 

0 

0 

0 

0 

0 

0 

1/G12 

The  compliance  matrix  as  shown  in  (7)  consists  of  eight  con¬ 
stants.  To  retain  the  five  independent  constants,  the  following 
relationships  are  required: 


v12  _  V21 

E1  "  E2  (8) 

v23  =  v32  (9) 

E2 

G23  =  2(1  +  v23)  00) 


Since  the  compliance  matrix  is  the  inverse  of  the  modulus  matrix, 
we  can  express  the  modulus  matrix  in  terms  of  the  engineering  constants 
[22,  23]: 


(1  V23^E1 

'll  "  “  : 

v23  "  2v12v21 

=  v12E2 
12  1  "  v23  “  2v12v21 


(ID 

(12) 


k  +  G 


23 


(13) 
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23  k  "  G23 

(14) 

44  =  G23 

(15) 

C66,"  G12 


In  Equation  (17),  k  is  the  bulk  modulus  for  plane  strain. 


(16) 

(17) 

(18) 


The  Micromechanics  Equations 


The  micromechanics  theory  to  be  used  is  based  upon  an  idealized 
geometry  which  leads  to  modified  rule-of -mixtures  Equations  [13].  The 
method  assumes  a  concentric  cylinder  model.  Figure  3  ,  in  which  the 
fiber  phase  is  embedded  in  matrix  material.  This  is  surrounded  by 
an  infinite  homogeneous  medium  indistinguishable  from  the  composite. 

The  fiber  is  considered  continuous  and  perfectly  aligned.  Both  fiber 
and  matrix  phases  are  assumed  homogeneous.  For  an  isotropic  matrix 
material  and  a  transversely  isotropic  fiber,  the  composite  engineering 
properties  are: 


(20) 
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where 


where 


where 


f  .  m 
v  +  n  v 

S'2  =  ^Z 

G12  Gm 


rm 

n  =  1/2(1  +  ^  ) 

<2 


Gm  = 


-m 


2(1  +  vm)  . 


23  f 


m 


v  nrv 

rr  +  -V 

G23  Gm 


=  3  -  4vm  +  G23 
^  4(1  -  vm) 


v  +  nkv 


m 


23 


m 


+  "k 


Hi,V 


'23 


m 


<m  = 


' 23 


-m 


2(1  -  2vm)(l  +  vm) 


■22 


2(1  -  v23  -  ^]2v2]  ) 


1  +  k-f 
k23 


2  (1  -  v”) 


=  4k23G23 
2  k23  +  yG23 


(21) 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


(30) 
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=  k23  "  yG23 
23  k23  +  hG23 


where  u  =  1  + 


2 

23  V1 2 


4k0~v 


'I 


(31) 

(32) 


and  the  coefficients  of  thermal  expansion  are 

«1  -  r  fefoV  +‘«W]- 

«2  =  (a2  +  “I  ^2^  +  (^  +  vm)amvm  -  v, po. 


(33) 


'12ul  (34) 

In  the  above  equations,  the  superscripts  f  and  m  are  fiber  and  matrix 
properties  and  v  is  the  volume  fraction. 


Constituent  Material  Properties 


As  presented  in  the  last  section,  the  laminate  will  be  composed  of 
transversely  isotropic  fibers  in  an  isotropic  matrix.  For  the  general¬ 
ized  plane  strain  finite  element  modeling  to  be  done,  a  full  set  of 
three-dimensional  engineering  properties  is  required. 

Since  this  investigation  was  planned  as  an  analytical  study,  the 
author  necessarily  relied  upon  existing  experimental  work.  Initially, 
the  extensive  experimental  characterization,  including  static,  visco¬ 
elastic,  and  fatigue  testing  at  various  evnironments  of  the  AS/3501-6 
graphite/epoxy  system  performed  by  Renton  and  Ho  [10,16,17]  seemed  ade¬ 
quate.  But  the  characterization  of  consituent  properties  was  limited 
in  two  areas.  First,  the  fibers  were  assumed  to  be  isotropic  and  sec¬ 
ond,  only  evaluation  of  the  in  plane  material  properties  were  consid¬ 
ered.  These  shortcomings  leave  this  characterization  inadequate  for 
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this  investigation. 

The  required  constituent  material  properties  to  be  considered  in 
this  investigation  are  presented  in  Table  1  and  were  obtained  from 
Reference  [20]. 

Laminate  properties  are  presented  in  Table  2.  Experimentally  de¬ 
termined  laminate  properties  were  obtained  from  Monib  and  Adams  [24] 
and  Renton  and  Ho  [16].  Monib  and  Adams  reported  =  0.88  x  10"6/°C. 
Upon  investigation,  this  value  was  found  to  be  unreasonably  high  and 

C 

was  replaced  with  =  0.15  x  10"  /°C  from  Reference  [25].  The  last 
column  of  Table  2  are  the  material  properties  determined  by  micro¬ 
mechanics  using  the  constituent  properties  from  Table  1. 

Due  to  the  very  low  predicted  values  for  the  major  Poisson's 
ratio,  v-j2,  obtained  by  the  micromechanics,  a  more  detailed  expres¬ 
sion  was  considered  to  improve  the  current  rule-of-mixtures  expres¬ 
sion.  From  Hashin  [22],  the  following  expression  for  v-j2  was 
considered: 


— 

* 

v12  = 

v12  + 

v12 

f  f 

,  m  m 

V1 2  = 

V12v 

+  v  v 

C 

1 

mw  1 

\  m 

* 

\  1  “ 

k23 

V1 2 

..IT 1 

.  f 

1  x  f  m 
T“)v  v 


23 


v  +  v_  +  1_ 
k23  k23  ^ 


(35) 

(36) 


(37) 


m 

where  k22,  k22  are  the  plane  strain  bulk  moduli  defined  in  equations 
(27)  and  (28).  In  Equation  (35)  the  term  v-|2  is  the  rule-of- 
mixtures  expression  previously  used.  The  second  term,  v*2,  is  an 


Table  1.  Constituent  Material  Properties  From 
Reference  [20],  for  Room  Temperature, 
Dry  Conditions. 


Hercules  AS- 
Graphite  Fiber 

Hercules  3501-6 
Epoxy  Matrix 

E-j  ( GPa) 

220. 

4.27 

E2(GPa) 

13.9 

4.27 

G12(GPa) 

34.5 

1.58 

G23(GPa) 

5.5 

1.58 

v12 

0.2 

0.34 

v23 

0.25 

0.34 

a-j  (10~6/°C) 

-0.36 

40.0 

a2(10'6/°C) 

18.0 

40.0 

Table 

2.  Comparison  of 

Laminate  Properties. 

Experimental,  Experimental, 

References  [24,25]  Reference  [16] 

Micromechanics 

Properties** 

(GPa) 

129.3 

130.7 

139.1 

E2(GPa) 

10.7 

6.19 

8.58 

G^2(GPa) 

6.65 

6.50 

5.84 

G23(GPa) 

3.54 

3.16 

v12 

0.31 

0.35 

0.25 

v23 

0.49 

0.36 

(10“6/°C) 

0.15* 

.11 

a2(10_6/°C) 

30.6 

31.23 

vf(%) 

62.0 

62.0 

*the  only  material  property  from  Reference  [25] 
**  using  constituent  properties  from  Table  1. 
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interaction  term  which  accounts  for  secondary  effects.  Note  that 

becomes  negative  when  vm  >  v[2  and  k£3  >  kjj3  which  is  the  case 

for  the  constituent  materials  in  Table  1;  in  fact,  vJ2  =  -.0037. 

* 

Since  v^2  is  small  and  negative,  the  rule-of-mixtures  expression 
originally  used  will  <be  considered  sufficient  to  predict  v^. 

At  this  point,  experimental  work  seems  necessary  to  resolve  the 
question  of  validity  of  the  micromechanics  and  the  constituent  material 
properties.  Since  the  present  study  is  primarily  concerned  with  the 
development  of  an  analytical  model ,  and  considerable  time  is  required  to 
perform  the  required  experimental  work,  the  micromechanics  shall  be 
assumed  valid.  The  material  properties  have  been  adjusted  to  those 
shown  in  Table  3.  A  new  comparison  of  the  predicted  material  proper¬ 
ties  from  micromechanics  to  experimental  values  is  shown  in  Table  4. 
The  predicted  material  properties  are  in  close  agreement  with  the  ex¬ 
perimental  values  of  References  [24,25].  The  selection  of  the 
constituent  material  properties  was  done  by  trial  and  error.  No 
specific  criterion  was  used  in  selecting  the  constituent  material 
properties  with  the  exception  of  G23  of  the  fiber.  The  transverse 
shear  modulus,  G23,  of  the  laminate  was  found  to  be  very  insensitive 
to  G23  of  the  fiber;  therefore,  G23  of  the  fiber  was  assumed  to  be 
60GPa,  an  amount  reported  by  Reference  [16]  as  an  upper  limit. 


Table  3.  Revised  Constituent  Material  Properties. 


Hercules  AS 
Graphite  Fiber 

Hercules  3501-6 
Epoxy  Matrix 

E-j  (GPa) 

206.0 

4.6 

E2(GPa) 

18.0 

4.6 

G]2(GPa) 

60.0 

1.65 

G23(GPa) 

6.77 

1.65 

v]2 

0.265 

0.39 

v23 

0.33 

0.39 

ai(io6/°c) 

0.39 

40.0 

a2(106/°C) 

20.0 

40.0 

Table  4.  Comparison  of  Laminate  Properties  With  Revised 
Constituent  Material  Properties. 


Experimental, 
References  [24,25} 

Experimental, 

Reference 

Micromechanics 

Properties* 

E-|  (GPa) 

129.3 

130.7 

129.5 

E2(GPa) 

10.7 

6.19 

10.4 

G12(GPa) 

6.65 

6.50 

6.35 

G23(GPa) 

3.54 

3.53 

V1 2 

0.31 

0.35 

0.31 

v23 

0.49 

0.47 

(106/°C) 

0.15** 

0.16 

a2(106/°C) 

30.6 

33.4 

vf(%) 

62.0 

62.0 

♦  using  constituent  properties  from  Table  3. 

**the  only  material  property  from  Reference  [25]. 


III.  STRAIN  ENERGY  RELEASE  RATE 


The  use  of  strain  energy  release  rate  as  a  fracture  growth  cri¬ 
terion  for  cracks  in  brittle  elastic  materials  is  one  of  long  standing. 
In  1921,  Griffith  [26]  noted  that  when  a  crack  was  introduced  to  a 
stressed  elastic  plate,  an  energy  balance  must  be  maintained  between 
the  decrease  in  the  elastic  strain  energy  and  the  increase  in  the 
surface  energy  resulting  from  growth  of  the  crack.  Griffith's  crack 
growth  criterion  for  fixed  grips  can  be  stated  as 


dU  _  dW 
da  da 


(38) 


where  U  is  the  elastic  strain  energy,  W  is  the  energy  required  for 
crack  growth,  and  a  is  the  crack  length.  The  term  is  usually  called 
the  strain  energy  release  rate,  G.  The  strain  energy  release  rate  is 
a  measure  of  the  driving  force  the  body  contains  to  extend  a  crack. 

The  term  ^  is  the  energy  consumed  during  crack  growth.  This  is  some¬ 
times  called  the  crack  resistance,  R.  The  units  associated  with  G  and 
R  are  energy  per  unit  of  crack  surface. 

Irwin  [27]  considered  an  infinite  plate  with  fixed  ends  containing 
a  crack  size  a.  If  forces  are  applied  to  the  crack  edge,  sufficient 
to  close  the  crack  over  an  infinitesimal  distance  (see  Figure  4)  work 
will  be  done.  If  this  process  is  reversed,  then  work  will  be  released 
from  the  plate.  This  work  per  unit  of  cracked  area  is  another  defini¬ 
tion  of  strain  energy  release  rate. 

Referring  once  again  to  Figure  4,  let  Au  be  a  vector  representing 
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the  relative  normal  and  sliding  displacements  between  the  crack  faces. 
Let  a  be  a  vector  defining  the  normal  and  shear  force  components  per 
unit  area  which  are  required  to  close  the  crack.  Then  Irwin's  crack 
closure  integral  can  be  expressed  as  : 


n  1  /s 

o  2Aa  /  # 


au  da 


By  writing  the  inner  product  in  component  form,  the  strain  energy 
release  rate  can  be  expressed  as  : 


where  : 


G  Gj  +  Gjj  +  Gin 


Gj  =  lim 

1  Aa  -*•( 


1  rL 
im  1  j 

->  O  2Aa  J 


aw  da 


/Aa  .  . 

Aw  Ha  (42) 

o„,  av  da 
yz 

o 

Gm  =  lim  1  /a  Ai)  .a  (43) 

111  Aa-o  IKZj  ^xz  Au  da 
o 

The  components  Gj,  Gjj,  and  Gj  j  j ,  also  referred  to  as  mode  I,  mode  II, 
and  mode  III  strain  energy  release  rates,  are  illustrated  in  Figure  5. 
These  correspond  to  the  crack  opening,  sliding,  and  tearing  modes, 
respectively. 


Finite  Element  Evaluation  of  Crack  Closure  Integral 

Rybicki  and  Kanninen  [5]  presented  a  very  simple  method  of  evalu¬ 
ating  the  Irwin's  crack  closure  integral  (39).  This  technique  greatly 
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simplifies  the  energy  release  rate  computation  because  knowledge  of 
the  singular  stress  field  near  the  crack  tip  is  not  required. 

Referring  to  Figure  6  (a),  first  the  nodal  forces  are  evaluated 
at  nodes  b  and  c  for  a  crack  of  size  a.  Then  nodes  b  and  c  are 
released  and  the  crack  is  allowed  to  extend  an  amount  Aa  as  in  Figure 
6  (b).  After  the  crack  extension  the  relative  displacements  between 
nodes  b1  and  e1  and  between  c'  and  f1  are  measured.  The  energy  release 
rate  can  now  be  expressed  as  the  work  required  to  close  the  crack  an 
amount  Aa  or  simply: 

GI  =  [Fzb(wb'  ’“e')  +  F2c(V  ‘wf,)]  (44) 

GII  =  2Za  ^Fyb(vb'  “V*  +  Fyc  "  vf'^  <45) 

GIII  =  2Za  ^Fxb^ub' "  ue'^  +  Fxc^uc' "  uf'^  (46) 

where,  for  example,  F  ^  is  the  nodal  force  in  the  z  direction  at  node 
b.  Observe  that  nodal  forces  at  node  b  are  composed  of  contributions 
from  both  elements  1  and  2.  One  expression  from  which  the  nodal  forces 
can  be  evaluated  for  the  linear  case  is: 

R®  =  q®  -Q®  1.  j  =  1,  NDOFE  (47) 

where  : 

0 

=  nodal  force  of  element 
K..  =  element  stiffness  matrix 

*  J 

qe.  =  element  nodal  displacements 

J 

6 

Q  .  =  applied  elemental  nodal  forces 
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NDOFE  =  number  of  degrees  of  freedom 
per  element. 


Delamination  Analysis  Based  Upon  Classical  Lamination  Theory 

In  this  section  a  method  of  delamination  analysis  will  be  presented 
based  upon  classical  lamination  theory  (CLT).  This  method  was  first 
used  by  O'Brien  [8]  for  the  case  of  a  constant  axial  strain  or 
"fixed-grip"  condition.  The  following  analysis  will  be  for  the  case  of 
constant  axial  stress. 

The  first  step  is  to  develop  the  rule  of  mixtures  expression  for 
the  change  in  the  axial  stiffness  which  is  common  to  both  load  cases. 
From  CLT  the  axial  stiffness  of  a  laminate.  Figure  7(a),  can  be 
expressed  as  : 


"LAM  a-|-|t 


where  a^  is  the  first  element  of  the  extensional  compliance  matrix 
and  t  is  the  total  laminate  thickness.  Now  consider  one  or  more 
planar  delaminations  to  extend  across  the  laminate.  These  delamina¬ 
tions  have  completely  severed  the  laminate  into  separate  sublaminates. 
The  sublaminates  are  assumed  to  be  gripped  in  such  a  manner  as  to  allow 
identical  axial  strainsin  all  the  sublaminates,  but  the  transverse 
strain  is  free  to  do  as  it  will.  The  axial  stiffness  after  the  full 
delamination.  Figure  7(b),  can  be  expressed  as: 


e*  =  i£  Ej  t, 

i=l 


(49) 


igure  7  .  Rule-of-mixtures  analysis  for  stiffness  loss. 
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where  =  the  axial  stiffness  of  the  ith  sublaminate  formed 
t.  =  the  thickness  of  the  ith  sublaminate  formed 
t  =  the  total  laminate  thickness 
n  =  the  number  of  sublaminates  formed. 

Since  the  axial  stiffness  of  the  sublaminate  is: 


then  E*  is  simply: 


11  Vi 


i=i  (an}i 


(50) 
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The  limiting  cases  for  the  axial  stiffness  have  now  been  established. 
These  are  the  totally  bonded  laminate  and  the  totally  delaminated 
laminate.  If  a  linear  variation  of  axial  stiffness  with  respect  to 
the  nondimensional  delamination  size  a/b,  as  in  Figure  7  (c),  is  assumed 
then  a  rule-of-mixtures  expression  is  derived  for  the  axial  stiffness, 

E  =  <e*-W!+elam  <52> 

Note  that  the  delamination  is  considered  to  extend  the  entire  length  of 
the  laminate  and  delaminations  are  considered  to  extend  inward  from  both 
sides  simultaneously. 

O'Brien  [8]  experimentally  verified  the  rule-of-mixtures 
expression  for  the  axial  stiffness  for  a  [+30/+30/90/90]s  laminate. 

But  the  verification  required  neglecting  the  contributions  of  the 
bending-extension  coupling  matrix.  Recall  that  the  axial  stiffness  of 
a  sublaminate  can  be  expressed  as: 


where  [a]  =  [A]'1  +  [A]''[B]  [D*]  [B]  [A]'1 
and  CD*]  =  [D]  -  [B]  W’CB] 

where  [A],  [B],  and  [D]  are  the  extensional,  coupling,  and  bending 
matrices  as  defined  in  reference  [28].  Provided  the  laminate  or  sub¬ 
laminate  in  question  is  symmetric,  a^  has  no  contribution  from  the  [B] 
matrix  since  it  vanishes  for  this  case.  But  for  the  laminate  examined 
by  O'Brien  delaminations  formed  at  the  -30/90  interfaces,  and  E*  was 
calculated  based  upon  the  following  three  sublaminates:  [±302]j,  [90^1T 
and  [f302]T.  Two  of  the  sublaminates:  [+302]  and  [+302]  are  nonsym- 
metric,  hence  the  coupling  matrix  [B]  is  not  zero.  O'Brien  neglected 
the  effect  of  the  [B]  matrix  to  obtain  the  best  match  to  experimental 
results.  The  reason  for  neglecting  the  [B]  matrix  is  due  to  the  fact 
that  instead  of  two  planar  delaminations  which  were  assumed  in  the 
model,  there  was  one  delamination  which  jumps  back  and  forth  from  one 
-30/90  interface  to  the  other  as  in  Figure  1(a).  The  delaminations 
were  believed  to  grow  in  this  manner  to  reduce  the  effect  of  the  bend¬ 
ing-extension  coupling  that  would  have  been  present  in  a  planar  dela¬ 
mination.  At  this  point,  it  cannot  be  concluded  that  the  effect  of 
the  [B]  matrix  in  general  can  be  neglected  without  further  investiga¬ 
tion  of  other  laminate  configurations. 

The  strain  energy  release  rate,  G,  for  a  planar  flaw  of  area.  A, 


can  be  expressed  as  [8], 
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where  W  is  the  potential  energy  of  the  applied  tractions  and  U  is  the 
elastic  strain  energy.  For  a  laminate  under  uniaxial  stress  a(in  which 
a  =  axial  force/area)  the  strain  energy  is: 


(62) 


34 


From  equations  (54),  (60),  and  (62)  we  obtain. 


g2t^ELAM  ~  E*) _ 

2[(e*-elam)  f-  +  Elam]2 


(63) 


Note  that  the  strain  energy  release  rate  becomes  infinite  as  E*  goes 
to  zero  and  goes  to  one. 
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IV.  APPLICATIONS 

In  this  chapter,  several  problems  will  be  considered  using  the 
analytical  tools  developed  in  the  previous  chapters.  An  extensive  dis¬ 
cussion  of  the  finite  element  method  used  in  this  chapter  is  contained 
in  Appendix  A.  First,  an  investigation  of  the  axial  laminate  stiffness 
as  a  function  of  a  reduction  in  the  matrix  modulus  is  conducted.  Once 
this  is  obtained,  the  matrix  degradation  as  a  function  of  fatigue  cycles 
is  established.  The  chapter  is  then  closed  with  a  parametric  study  of 
the  effect  of  matrix  degradation  on  the  available  strain  energy  release 
rate  for  a  delamination  at  the  -45/90  ply  interface.  Throughout  this 
entire  chapter,  only  a  [±45/902^  laminate  will  be  analyzed.  Residual 
stresses  induced  by  moisture  and  thermal  effects  have  been  neglected. 

Laminate  Stiffness  as  a  Function  of  Matrix  Modulus 

In  this  section,  the  axial  laminate  stiffness  as  a  function  of 
matrix  modulus  for  a  [±45/9023s  test  specimen  will  be  investigated. 

The  main  point  of  this  will  be  to  assess  the  accuracy  of  classical 
lamination  theory  (CLT)  versus  a  finite  element  method  (FEM)  which 
accounts  for  free  edge  effects.  For  this  study,  the  micromechanics 
theory  previously  outlined  was  used  with  the  constituent  properties 
from  Table  3.  The  test  specimen  geometry  used  for  the  laminate  was 
from  Reference  [10],  this  information  being  contained  in  Appendix  B. 

Recall  from  the  previous  section,  from  CLT  the  axial  stiffness 
for  the  laminate  can  be  expressed  as: 

lLAM  ant 


(64) 
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where  is  the  first  term  in  the  extensional  compliance  matrix  and 
t  is  the  laminate  thickness.  For  the  FEM  solution,  the  axial  laminate 
stiffness  can  be  calculated  as: 


(65) 


where  om  =  Jaxx  dydz.  (66) 

The  term  can  easily  be  identified  as  the  average  axial  stress 

and  e  is  the  axial  strain. 

o 

Since  the  free  edge  stresses  are  local  to  the  free  edge,  only  a 
portion  of  the  entire  width  of  the  laminate  was  modelled  using  the 
finite  element  method.  The  width  of  the  laminate  selected  to  be  model¬ 
led  was  determined  by  the  distance  required  from  the  free  edge  for 
the  stresses  to  converge  to  the  CLT  solution.  Hence,  the  elasticity 
solution  was  obtained  by  a  combination  of  the  FEM  solution  near  the 
free  edge  and  the  CLT  solution  for  the  interior  of  the  laminate. 

Because  of  this,  the  expression  for  the  axial  laminate  stiffness  had  to 
be  modified.  It  is  a  simple  matter  to  show  that  the  axial  stiffness 
for  the  elasticity  solution  is  a  rule-of-mixtures  expression  obtained 
from  the  two  solutions,  if  one  considers  two  parallel  bars  which  are 
deformed  uniformly.  The  rule-of-mixtures  expression  can  be  stated  as: 

ELAM  =  ECLT  ACLT  +  EFEM  afem  ^ 

where  E^y  and  Ep^  are  the  axial  stiffnesses  from  the  CLT  and  the  FEM 
solutions,  respectively,  and  A^y  and  ApEM  are  the  fractions  of  the 
cross-sectional  area  for  the  respective  solutions.  The  elasticity 
solution  is  depicted  in  Figure  9.  Due  to  symmetry,  only  one-quarter 


(a)  Cross-section  of  laminate. 
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Region  of  classical 
lamination  solution. 


Region  of  free  edge 
stress,  FEM  solution. 


(b)  One-quarter  of  cross-section  analyzed  due 
to  symmetry  conditions. 
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(c)  Only  free  edge  region  requires  finite  element  modelling. 


Figure  9.  Schematic  of  elasticity  solution. 
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of  the  laminate  required  modelling  and  of  that  quarter,  only  the  portion 
to  the  right  of  the  z'axis  required  finite  element  modeling. 

The  axial  stiffness  was  evaluated  for  six  values  of  the  matrix 
modulus,  from  the  initial  matrix  modulus  of  4.6  GPa,  down  to  a  value 
of  0.0723  GPa.  This  wide  range  of  values  for  the  matrix  modulus 
required  using  two  different  meshes  for  the  finite  element  model.  The 
first  mesh,  see  Figure  10,  was  selected  to  model  one-third  of  the 
laminate  half-width  and  maintain  a  maximum  element  aspect  ratio  of  10:1. 
This  mesh  was  adequate  for  the  initial  matrix  modulus  and  63%  of  the 
initial  matrix  modulus,  but  for  25%  of  the  initial  matrix  modulus,  the 
stresses  along  the  z'  axis  did  not  converge  to  the  CLT  solution;  hence, 
the  region  of  free  edge  effects  was  not  modelled  sufficiently.  To 
correct  this,  two-thirds  of  the  laminate  half-width  was  modelled.  Due 
to  program  limitations,  additional  elements  were  not  available,  so  the 
first  mesh  was  simply  expanded.  The  second  mesh  is  shown  in  Figure  11  . 

In  expanding  the  mesh,  the  aspect  ratio  increased  proportionately  to 
20:1.  The  second  mesh  was  used  to  compute  the  axial  stiffness  for 
matrix  modulus  values  from  25%  to  1.6%  of  the  initial  matrix  modulus. 

Both  meshes  contained  120  elements  and  1158  degrees  of  freedom. 

Presented  in  Figure  12  is  the  axial  laminate  stiffness  as  a  func¬ 
tion  of  matrix  modulus.  In  evaluating  the  laminate  stiffness,  only  the 
matrix  modulus  was  allowed  to  vary  while  all  other  constituent  material 
properties  were  held  constant.  On  the  x-axis,  the  matrix  modulus  was 
normalized  by  ti al  *  w^cl1  was  the  initial  matrix  modulus  from 
Table  3.  On  the  y-axis,  the  laminate  stiffness  was  normalized  by 


E 


REF5 


which  was  defined  as  the  laminate  stiffness  obtained  from  the 


1  ply  thickness 


1  ply  thickness 


function  of  matrix  modulus. 
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elasticity  solution  when  Em/Elj1nl.t1.al  =  1.0.  The  value  of  EREF  was 
found  to  be  23.02  GPa.  It  should  be  noted  that  EREF  will  be  used 
throughout  the  rest  of  this  work  in  presenting  the  results  in  a  nondi¬ 
mens  ional  form.  From  Figure  12,  the  CLT  solution  is  shown  to  be  stif- 
fer  than  the  FEM  solution  and  for  smaller  values  of  the  matrix  modulus 
the  discrepancy  becomes  quite  considerable.  This  discrepancy  is  due 
to  the  fact  that  as  the  matrix  modulus  goes  to  zero,  CLT  yields  a 
finite  value  of  12.7  GPa  for  the  laminate  stiffness  while  the  finite 
element  stiffness  appears  to  vanish. 

Presented  in  Figure  13  is  the  error  in  the  laminate  stiffness  as  a 
function  of  matrix  modulus  when  calculated  by  CLT.  The  x-axis  is  the 
normalized  matrix  modulus  as  used  previously.  For  the  y-axis,  E^y 
and  E^y^yare  the  axial  laminate  stiffnesses  obtained  by  CLT  and  the 
elasticity  solution,  respectively.  It  can  be  seen  for  the  initial 
value  of  the  matrix  modulus  that  the  difference  is  a  modest  2.4%,  but 
for  the  normalized  values  of  the  matrix  modulus  of  0.1  and  0.016,  the 
difference  increased  to  14.7%  and  48.0%,  respectively. 

The  reason  for  this  extreme  discrepancy  in  the  CLT  solution  is 
that  the  shear  modulus  of  the  laminate  tends  to  maintain  a  value 
on  the  same  order  of  magnitude  as  the  matrix  modulus.  Hence,  as 
the  matrix  modulus  goes  to  zero,  so  does  the  shear  moduli.  This  result¬ 
ing  decrease  in  the  shear  modulus  allows  for  considerable  shear  defor¬ 
mation  which  tends  to  warp  the  laminate  cross-section.  Due  to  the  fun¬ 
damental  Kirchhoff-Love  hypothesis  of  CLT,  it  is  unable  to  model  the 
necessary  shear  deformation  encountered  with  the  very  small  shear  modu¬ 
lus,  and  this  accounts  for  the  erroneously  stiff  solutions  predicted  by 
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CLT.  It  should  be  noted  that  the  FEM  used  in  analyzing  the  laminate 
contained  sufficient  freedom  to  model  the  shear  warping. 

Matrix  Degradation  as  a  Function  of  Fatigue  Cycles 

In  this  section,  we  shall  evaluate  the  matrix  modulus  as  a  function 
of  the  fatigue  cycles.  This  will  be  based  upon  the  mean  strain  fatigue 
test  data. from  Ho  and  Schapery  [10]  which  is  included  in  Appendix  B. 
Throughout  this  section,  the  laminate  will  be  considered  without  any 
delaminations.  The  study  of  the  effect  of  delaminations  will  be  con¬ 
ducted  in  the  next  section. 

Constant  stress  amplitude  fatigue  tests  were  conducted  with  maxi¬ 
mum  tensile  stresses  of  50%,  55%,  and  60%  of  the  ultimate  tensile 
stress  of  the  laminate.  The  ratio  of  the  minimum  stress  to  maximum 
stress  was  +0.1.  Additional  information  concerning  test  conditions 
are  contained  in  Appendix  B. 

Presented  in  Figure  14  is  the  mean  strain  response  as  a  function 
of  the  number  of  fatigue  cycles  for  each  of  the  tests.  These  curves 
are  somewhat  idealized  due  to  scatter  and  the  limited  amount  of  test 
data.  A  further  discussion  of  the  mean  strain  idealization  is  con¬ 
tained  in  Appendix  B.  The  mean  strain  behavior  is  characterized  by 
an  increasing,  predominately  linear,  response  in  logarithmic  scales. 

Near  failure,  the  mean  strain  breaks  away  from  its  previous  linear 
path  and  increases  very  rapidly  to  failure. 

Throughout  the  remainder  of  this  work  the  fatigue  tests  shall  be 
analyzed  as  though  they  were  static  tests.  Simply  stated,  no  dynamic 
effects  will  be  considered.  Also,  the  effects  of  the  stress  amplitude 


t 


Figure  14.  Idealized  mean  strain  cu 
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will  not  be  addressed  '/directly.  Instead,  only  the  mean  stress  input 
and  the  mean  strain  response  will  be  addressed.  From  this  very  nar¬ 
row  viewpoint,  the  fatigue  tests  somewhat  resemble  creep  tests.  With 
this  in  mind,  we  shall  define  the  apparent  laminate  stiffness  as: 

E  APPARENT  =  -  (68) 

e 

where  a  is  the  mean  stress  and  F  is  the  mean  strain.  The  apparent 
laminate  stiffness  for  the  three  fatigue  tests  is  shown  in  Figure  15. 
The  apparent  stiffness  was  normalized  by  the  term  EREF,  which  is  the 
initial  laminate  stiffness  previously  defined.  From  this  figure  it  is 
easy  to  see  that  in  all  of  the  tests,  the  stiffnesses  degraded  to  ap¬ 
proximately  75%  of  their  original  stiffness  by  the  tenth  cycle,  which 
is  quite  interesting.  One  possible  explanation  for  the  30%  reduction 
in  the  axial  stiffness  is  the  initial  matrix  damage  as  a  result 
of  the  residual  thermal  stresses  incurred  during  the  cool  down  from 
the  laminate  cure  temperature. 

The  relationship  of  matrix  modulus  to  fatigue  cycles  was  difficult 
to  obtain  due  to  the  lack  of  explicit  functions.  This  resulted  in 
determining  the  desired  function  graphically.  First,  the  elasticity 
solution  for  the  laminate  stiffness  as  a  function  of  matrix  modulus 
from  Figure  12  was  inverted  to  yield  the  matrix  modulus  as  a  function 
of  laminate  stiffness.  Then,  the  apparent  laminate  stiffness  as  a 
function  of  the  fatigue  cycles  from  Figure  15  were  substituted  in  for 
the  laminate  stiffness  in  the  previous  function  which  then  yielded 
the  matrix  modulus  as  a  function  of  the  fatigue  cycles.  The  procedure 
is  depicted  in  Figure  16.  The  matrix  modulus  as  a  function  of  the 
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Figure  15.  The  apparent  laminate  stiffness  as  a  function 
of  fatigue  cycles. 
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Inverting  (a)  yields. 


(b)  Em  =  Em(EUm) 
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Substituting  (c)  into  (b),  with  EAppARENT  =  ELAM,  yields. 


Figure  16.  Analytical  versus  graphical  method  for  obtaining 
matrix  modulus  as  a  function  of  fatigue  cycles. 
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number  of  fatigue  cycles  is  presented  in  Figure  17.  Since  the  curves 
are  represented  by  only  four  points,  considerable  detail  is  lost.  This 
is  particularly  unfortunate,  since  the  data  points  at  the  matrix  modulus, 
ratio  values  of  .25  and  .1  bracket  the  location  where  the  mean  strain 
began  its  sharp  upward  turn  which  would  appear  to  be  a  critical  point 
of  concern.  For  this  reason,  matrix  modulus  values  of  .25  and  .1,  in 
addition  to  1.0,  will  be  considered  for  a  parametric  study  of  interply 
delamination  in  the  next  section. 

i 

Strain  Energy  Release  Rate  as  a  Function  of  Matrix  Degradation 

In  this  section,  a  parametric  study  of  the  effect  of  matrix  soften¬ 
ing  upon  the  strain  energy  release  rate  as  a  function  of  delamination 
length  is  made.  Only  three  values  of  the  matrix  modulus  are  studied» 
because  of  the  considerable  amount  of  computational  time  involved  in 
the  finite  element  analysis.  The  three  values  of  matrix  modulus  ratio, 
^/Einitial  ’  ^died  are  1»  0*25,  and  0.1.  Matrix  modulus  ratios  of 
0.25  and  0.1  bracket  the  location  of  the  breakaway  point  on  the  mean 
strain  versus  fatigue  cycle  curves.  Further  discussion  of  this  point 
is  contained  in  the  previous  section. 

The  selection  of  the  most  energetically  favorable  location  for  the 
free  edge  delamination  was  done  using  O'Brien's  method  of  analysis. 

The  analysis  indicated  that  more  energy  was  available  for  delamination 
at  the  -45/90  interface  than  at  any  other  point.  This  conclusion  is 
supported  by  the  study  of  Law  [29].  In  a  similar  laminate,  a  [±25, 

^2-^s  stackln9  sequence.  Law  used  the  finite  element  method  to  predict 
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Figure  17.  Matrix  modulus  as  a  function  of  fatigue  cycles. 
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delamination  at  the  -25/90  interface.  This  interface  corresponds  to 
the  -45/90  interface  of  the  laminate  considered  in  this  study.  Hence, 
all  the  subsequent  delamination  analysis  was  conducted  at  the  45/90  ply 
interface;  examination  of  [±45/902]$  test  specimens  seems  to  support 
this  assumption. 

The  finite  element  modeling  was  performed  with  the  same  meshes  as 
were  used  to  evaluate  the  axial  stiffness.  For  the  matrix  modulus 
ratio  of  one,  the  strain  energy  release  rate  was  calculated  for  a 
series  of  nine  crack  extensions  for  both  meshes.  The  strain  energy 
release  rate,  calculated  by  the  crack  closure  technique,  displayed  con¬ 
siderable  insensitivity  to  the  mesh  coarseness.  The  difference  in  the 
coarse  mesh  prediction  compared  to  the  fine  mesh  varied  from  5%  at  one 
ply  thickness  from  the  free  edge  to  2%  at  four  ply  thicknesses.  It 
was  decided  to  use  the  coarse  mecn  of  Figure  11  for  all  the  analysis. 
For  the  same  number  of  crack  extensions,  the  energy  release  rate  could 
be  evaluated  for  a  length  of  10  ply  thicknesses  from  the  free  edge 
versus  only  4  ply  thicknesses  with  the  finer  mesh  at  a  sacrifice  of 
about  5%  error  in  the  accuracy  near  the  free  edge. 

The  first  load  case  was  a  fixed  grip  condition.  The  energy  release 


rates  as  a  function  of  the  delamination  for  each  of  the  values  of 
Em/E^  ...  ,  are  presented  in  Figures  18,  19,  and  20.  The  length  of 

im  U  I  cl  I 

delamination,  a,  is  normalized  by  the  ply  thickness,  h,  throughout  the 


remainder  of  this  study.  The  strain  energy  release  rate,  G,  is  normal¬ 
ized  by  the  maximum  value  of  the  total  strain  release  rate  for  a  parti¬ 
cular  value  of  the  matrix  modulus  ratio.  These  figures  were  included 


to  illustrate  the  relative  contributions  of  the  individual  components 
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Figure  19.  Components  of  strain  energy  release 
rate  for  =  0.25. 
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of  the  strain  energy  release  rate.  The  mode  III  strain  energy,  Gjjj, 
was  typically  one  to  two  orders  or  magnitude  less  than  the  other  strain 
energy  release  rate  components  and  has  been  neglected.  The  mode  II 

i 

energy  release  ra*te  is  observed  to  decrease  as  the  matrix  modulus 
decreases.  This  trend  is  due  to  the  corresponding  decrease  in  the 
shear  moduli  which  restricts  the  ability  of  the  laminate  to  store 
energy  in  this  mode.  Hence,  less  energy  is  available  to  be  released 

in  mode  II.  This  is  summarized  by  Figure  21. 

The  total  strain  energy  release  rate,  G,  as  a  function  of  the 

normalized  delamination  length  is  shown  in  Figure  22.  The  total 
strain  energy  release  rate  was  normalized  by  the  applied  strain,  eQ, 
the  total  laminate  thickness,  t,  and  the  reference  laminate  stiffness, 
Eref*  The  strain  energy  release  rate  obtained  by  the  finite  element 
method,  shown  by  the  solid  curves,  appears  to  converge  to  the  values 
obtained  from  O'Brien's  method,  shown  by  the  dashed  lines. 

The  expression  for  the  energy  release  rate  under  a  fixed  grip  con¬ 
dition  for  O'Brien's  method  can  be  easily  obtained  from  equation  (63). 
Recall  that  the  expression  in  equation  (63)  is  for  G  under  a  specified 
stress  loading  condition  and  that  the  term  in  the  denominator  is 
the  laminate  stiffness;  viz: 

a  2t(ELAM-E*) 


2  2  2 
Simply  substituting  eQ  into  the  above  equation  for  a  /E  will  yield 

the  fixed  grip  expression  for  G,  hence: 


G  -  ”V  Ielam-e*> 


(70) 


21.  Mode  II  fraction  of  total  energy  release 
rate  as  a  function  of  crack  length. 
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From  the  previous  expression,  it  can  be  seen  that  O'Brien's  method  for 
the  fixed  grip  condition  is  independent  of  the  delamination  length. 

What  is  quite  remarkable  is  the  accuracy  of  this  expression  for  values 
of  the  matrix  modulus  ratio  of  0.1,  whereas  classical  lamination  theory 
predicted  a  14.7%  error  in  the  laminate  stiffness.  This  apparent 
discrepancy  can  be  accounted  for  because  the  energy  release  rate  is 
actually  proportional  to  the  term  E^^-E*  which  is  related  to  dE/da. 
This  seems  to  imply  that  the  linear  variation  of  the  laminate  stiffness 
as  a  function  of  delamination'  size  is  valid,  and  as  shown  by  Figure  23, 
the  lamination  theory  can  accurately  predict  dE/da.  For  the  fixed 
grip  loading  case,  this  will  allow  the  use  of  the  finite  element  method 
to  evaluate  the  strain  energy  release  rate  in  the  vicinity  of  the  free 
edge  and  O'Brien’s  lamination  theory  to  evaluate  the  strain  energy  in 
the  interior  of  the  laminate. 

Now  consider  the  implications  of  Figure  22  with  respect  to  the 
fracture  behavior  of  the  laminate.  Since  a  fixed  grip  (specified 
strain)  was  applied,  this  figure  could  represent  a  constant  strain 
amplitude  fatigue  test.  For  this  type  of  a  fatigue  test,  the  mean 
strain  input  would  be  constant  and  mean  stress  response  would  be  de¬ 
creasing.  The  mean  stress-strain  behavior  would  somewhat  resemble  a 
viscoelastic  relaxation  test.  At  this  point  we  shall  assume  that  the 
energy  required  for  delamination  is  constant,  and  shall  be  denoted  as 


G  ,  the  critical  energy  release  rate, 
c 


Referring  to  Figure  22,  at  the  beginning  of  the  fatigue  test  the 
available  energy  for  delamination  is  shown  by  the  curve  corresponding 


to  Em/Erj'nitiai  =  1.0.  As  fatiguing  continues,  the  matrix  becomes 


O'Brien's  method 
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damaged,  resulting  in  a  decrease  in  the  apparent  matrix  modulus.  In 
the  figure,  it  should  be  noted  that  as  the  matrix  modulus  decreases, 
the  available  strain  energy  release  rate  also  decreases.  Hence,  a 
critical  period  for  constant  strain  amplitude  fatigue  would  be  early 
in  the  test  when  the  available  strain  energy? release  rate  is  largest. 

Since  the* strain  energy  release  rate  as  a  function  of  delamina- 

* 

tion  size  is  an  increasing  function  and  does  not  exhibit  any  minimum 
points,  a  de lamination  will  grow  when  the  Griffith  criterion,  G  =  Gc, 
is  satisfied.  Referring  to  Figure  24,  the  intersection  point  of  the 
critical  strain  energy  release  rate,  G  ,  and  the  available  strain 

V 

energy  curve  shall  be  defined  as  ac»  the  critical  delamination  size. 
During  the  initial  phase  of  the  fatigue  test,  an  initial  flaw  in  the 
laminate  which  was  smaller  than  ac  would  remain  stable  and  not  grow. 
Initial  flaws  which  are  equal  to  or  larger  than  a  would  grow  due  to 
the  excess  energy  available  and  since  there  appears  to  be  no  minimum 
point  in  the  available  energy  the  growth  will  be  unstable.  Once 
unstable  delamination  has  begun,  the  delamination  will  cover  the  entire 
specimen  width  and  specimen  failure  would  likely  result. 

A  constant  stress  loading  condition  is  considered  next.  Once 
again,  the  finite  element  solution  and  the  classical  lamination  theory 
based  method  of  O'Brien  is  compared.  O'Brien's  method  has  been  pre¬ 
sented  in  a  previous  section.  At  this  point,  a  description  of  how  the 
fixed  grip  solution  obtained  from  the  finite  element  method  was  con¬ 
verted  to  a  constant  stress  condition  will  be  discussed. 

Since  linear  elastic  fracture  mechanics  has  been  employed,  the 
strain  energy  release  rate  is  proportional  to  the  square  of  the  applied 


Strain  energy  release  rate. 
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Figure  24.  Delamination  growth  criterion. 


62 


strain.  This  factor  is  employed  to  scale  the  constant  strain  values 
of  the  strain  energy  release  rate.  When  the  laminate  is  loaded  ini¬ 
tially  the  delamination  size  is  zero  and  the  stress  can  be  expressed 
as: 

%  =  initial  eo  (71) 

When  a  delamination  of  size  a  is  introduced,  the  laminate  stiffness 
decreases  to  a  value  of  E,  and  to  maintain  a  constant  stress  of  aQ,  the 
strain  must  be  increased  such  that: 


c  =  !mniai  (72) 

E  0 

In  the  above  expression,  the  term  E^^.  -]/E  is  the  amplification  fac¬ 
tor  required  to  maintain  a  constant  stress  condition.  Since  the  lami¬ 
nate  stiffness  is  a  function  of  the  delamination  size,  the  correction 
factor  is  best  defined  as  the  laminate  stiffness  evaluated  at  a  delami¬ 
nation  size  of  zero  divided  by  the  laminate  stiffness  evaluated  at  the 
current  delamination  size.  Hence,  to  obtain  the  strain  energy  release 


rate  for  a  constant  aQ,  simply  multiply  the  strain  energy  release 
rate  for  a  constant  eQ  by  the  square  of  the  amplification  factor  de¬ 
fined  above. 


Presented  in  Figure  25  is  the  total  strain  energy  release  rate 
as  a  function  of  the  normalized  delamination  length.  The  strain 
energy  release  rate  was  normalized  by  the  applied  stress,  aQ,  the 
total  laminate  thickness,  and  the  reference  laminate  stiffness,  ER£p 
The  strain  energy  release  rate  is  observed  to  increase  as  the  matrix 
modulus  decreases.  For  the  matrix  modulus  ratio  of  1,  there  is 
reasonable  agreement  between  O'Brien's  method  and  the  finite  element 
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Figure  25.  Total  normalized  strain  energy  release 
rate  as  a  function  of  delamination 
length/ply  thickness  for  a  constant 
stress  condition. 
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method,  but  this  deteriorates  as  the  matrix  modulus  decreases.  The 
discrepancy  of  O'Brien's  method  is  due  to  the  inability  of  the  lamina¬ 
tion  theory  to  predict  the  laminate  stiffness.  Despite  this  problem, 
O'Brien's  method  is  of  considerable  importance, because  it  contains  all 
of  the  essential  elements  required  to  describe  the  energy  release 
rate, and  substantial  insight  can  be  obtained  by  the  study  of  it. 

Consider  a  laminate  with  an  initial  flaw  size  of  aQ  and  a  critical 
strain  energy  release  rate  of  G  .  As  this  laminate  is  fatigued  in  a 
constant  maximum  stress  fatigue  test,  the  matrix  modulus  degrades, 
and,  as  a  result,  the  available  energy  increases.  This  process  is 
depicted  by  points  A,  B,  and  C  in  Figure  26.  Point  C  is  a  critical 
point  since  the  energy  available  is  equal  to  the  critical  strain  energy 
release  rate.  At  this  point,  the  initial  flaw  will  propagate  unstably 
across  the  entire  sample  width  since  the  available  energy  appears  to 
continue  to  increase.  Once  the  laminate  has  fully  delaminated, 
failure  will  occur  very  quickly  due  to  the  degraded  matrix  properties. 

Now  if  a  critical  energy  release  rate  is  assumed,  we  can  deter- 
mine  the  critical  flaw  size  from  Figure  25.  Assuming  Gc  =  175  N*M/M  , 
the  corresponding  flaw  sizes  for  the  50%,  55%,  and  60%  UTS  fatigue 
tests  are  2.5h,  1.5h,  and  Ih  respectively,  where  h  is  the  ply  thick¬ 
ness.  This  range  of  flaw  sizes  is  reasonable  for  an  edge  flaw  in  the 
type  of  composite  considered. 
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VI.  CONCLUSIONS 

Matrix  degradation  modelled  by  a  reduced  matrix  modulus  through 
mieromechanics  principals  shows  considerable  promise  as  a  means  to 
analyze  fatigue  behavior.  The  study  conducted  was  of  a  preliminary 
nature  in  that  only  mean  stress  and  mean  strain  fatigue  responses 
were  considered,  and  the  influence  of  stress  gradients  on  matrix 
damage  was  neglected.  The  major  purpose  of  this  study  was  to  assemble 
the  necessary  analytical  tools  required  to  analyze  the  combined  effects 
of  matrix  degradation  and  delamination.  Considerable  energy  was 
expended  in  developing  a  generalized  plane  strain  finite  element 
program  capable  of  analyzing  free  edge  effects  and  performing  fracture 
analysis.  Throughout  this  study,  the  accuracy  of  the  finite  element 
method  and  classical  lamination  theory  was  assessed. 

In  evaluating  laminate  stiffness,  it  was  found  that  as  the  matrix 
modulus  was  reduced,  considerable  shear  deformations  occurred  which 
lamination  theory  was  unable  to  model.  The  shear  deformations,  which 
are  associated  with  free  edge  effects,  were  found  to  penetrate  a  con¬ 
siderable  distance  in  from  the  free  edge.  For  the  smallest  value  of 
the  matrix  modulus  considered,  the  finite  element  results  showed  that 
significant  free  edge  effects  extended  five  laminate  thicknesses  into 
the  interior  of  a  [±4.5/9023s  laminate. 

A  parametric  study  of  free  edge  delaminations  as  a  function  of 
matrix  damage  was  performed  for  constant  maximum  strain  and  constant 
maximum  stress  fatigue  tests.  For  the  constant  maximum  strain 
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fatigue  test  the  available  energy  for  del  ami nation  was  greatest  at  the 
beginning  of  the  test  and  then  was  subsequently  reduced  by  the  matrix 
damage.  For  the  constant  maximum  stress  fatigue  test,  the  effect 
of  matrix  damage  was  to  increase  the  energy  available  for  delamination. 
Once  sufficient  energy  became  available,  a  flaw  was  shown  to  propagate 
unstably  across  the  entire  laminate  width. 

The  accuracy  of  the  classical  lamination  theory  based  delamination 
analysis  of  O'Brien  suffered  from  the  inaccurate  laminate  stiffness 
predictions  due  to  the  influence  of  the  free  edge.  This  problem  was 
amplified  due  to  the  small  width  of  the  laminate.  None  the  less,  con¬ 
siderable  insight  can  be  obtained  by  the  use  of  O'Brien  s  method, since 
it  contains  all  of  the  essential  elements  which  are  required  to 
evaluate  the  strain  energy  release  rate,  in  a  simple  closed  form  solu¬ 
tion.  The  only  major  limitation  of  O'Brien's  method  is  it  does  not 
provide  values  for  the  separate  modes  of  strain  energy  release  rate, 

Oj ,  Gjj,  Gjjj* 
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VII.  RECOMMENDATIONS  FOR  FUTURE  WORK 

In  order  to  develop  a  more  realistic  fatigue  damage  model  for 
fibrous  composites,  several  additional  parameters  must  be  included 
other  than  those  considered  in  this  study.  Influence  of  such  factors 
as  the  effect  of  stress  gradients,  the  stress  amplitude,  cyclic  fre¬ 
quency,  residual  stresses,  and  environmental  effects  all  deserve  atten¬ 
tion.  The  most  important  consideration  in  developing  an  analytical 
model  is  maintaining  a  firm  physical  foundation. 

At  this  point,  considerable  experimental  work  is  needed  to  under¬ 
stand  fatigue  behavior  and  the  sequence  of  events  which  lead  to  fail¬ 
ure.  Such  a  study  should  try  to  assess  the  relative  influence  of 
matrix  degradation  and  delamination.  Experimental  techniques  such  as 
dye  penetrant  enhanced  X-radiography  could  be  used  to  study  fatigue 
delamination  growth.  The  effects  of  stress  gradients  could  be  evalu¬ 
ated  by  conducting  tests  of  specimen  widths  in  an  effort  to  isolate 
the  free  edge  effects.  It  is  believed  that  until  a  sufficient  experi¬ 
mental  data  base  has  been  established,  it  will  be  difficult  to  estab¬ 
lished  an  accurate  fatigue  damage  model. 

One  last  remark  to  future  researchers,  "I  think  you  will  find  that 
if  you  plot  it  on  log-log  paper  it  will  be  linear,"  from  Reference 
[14]. 
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FINITE  ELEMENT  FORMULATION 

Consider  the  linear  elastic  lamina  shown  in  Figure  A1 .  Each 
lamina  is  discretely  modeled  as  a  homogeneous  orthotropic  material. 
The  modulus  matrix  in  the  principal  material  coordinates  (1,  2,  3), 
[21],  is 


The  modulus  matrix  will  be  expressed  in  terms  of  engineering 
constants  referenced  to  the  principal  material  coordinates.  The  nota¬ 
tion  and  definition  of  the  following  engineering  constants  was  adopted 
from  Jones  [28]: 

E, ,  E2,  E3  =  Young's  moduli  in  1,  2,  and  3  directions 
respectively. 

v..  =  Poisson's  ratio  for  transverse  strain  in  the 

1J  j-direction  when  stressed  in  the  i -direction. 

G--, ,  G, ?  =  Shear  moduli  for  2-3,  3-1,  and  1-2  planes 
J  J  respectively. 

The  modulus  matrix  can  be  expressed  in  terms  of  the  engineer¬ 


ing  constants  as: 
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Figure  Al.  Principal  material  coordinates,  (1,2,3)  used 
to  define  the  engineering  constants. 
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When  the  lamina  is  rotated  an  arbitrary  amount  about  the  3,  z  axis. 
Figure  A2,  the  modulus  matrix  is  of  a  monoclinic  material  form 
with  respect  to  the  global  x,  y,  z  coordinate  system.  The  monoclinic 
modulus  matrix  is 


(A-2) 


The  transformation  from  principal  material  coordinates 
to  the  global  x,  y,  z  coordinates  by  a  rotation,  e,  about  the  z  axis 
are  [21]: 


CH  =  m4Cn  +  2m2n2(C12  +  2Cgg)  +  n4C22 
C]2  =  m2n2(^ir  +  C22  -  4Cgg)  +  (m4  +  n4)C12 
C13  =  m  Ci3  +  n  C23 

C-jg  -  mn [m  C-j-j  -  n  C22  -  (m  -  n  )(C-j2  +  2Cgg)] 
C22  =  n%i  +  2m2n2(C12  +  2C66)  +  m4C22 
C23  =  n  C13  +  m  C23 

C2g  =  mn[n2C]1  -  m2C22  +  (m2  -  n2)(C]2  +  2Cgg)] 
C33  =  C33 

C36  =  mn^ci3  '  c23^ 

C44  =  m%4  +  n%5 
G45  =  mn^G55  " 

C55  =  m  C55  +  n^44 
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=  m2n2(C11  +  C22  =  2C-j 2 )  +  (m2  -  n2)2  Ce 


(A— 3 ) 


where  m  =  cos  e,  n  =  sin  e. 

Note  that  the  above  transformation  is  based  upon  the  right  hand 
rule  for  positive  rotations  about  the  z  axis,  Figure  A2.  Note  that  an 
opposite  convention  was  used  in  [22], 

The  following  displacement  field  is  assumed  [n], 
u(x,  y,  z)  =  e0x  +  U(y,  z) 
v(x,  y,  z)  =  V(y,  z) 
w(x ,  y,  z)  =  W(y,  z) 

Let  us  determine  the  class  of  problems  that  can  be  solved.  Differen¬ 
tiating  the  assumed  displacments  yields  the  following  strains  and 
curvatures. 
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r  t/2 

V  Mxy  %/2  (V  V  Txy>  zdz 

rt/2  2 

B,.,  D..  =/  Q,,(l,  z,  z‘)  dz 

J  J  -t/2  J 


=  s12  s16 

S12  S22  S26 
_S16  S26  S66 


(A— 5) 


Substituting  in  the  strains  and  curvatures  from  equation  (A-4)  into 
equation  (A-5)  we  obtain 
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The  most  important  implication  of  the  assumed  displacement  field  is 

that  the  normal  strain  ex  =  eQ  and  the  kx  and  kxy  curvatures  are  zero. 

Notice  the  effect  of  the  bending-extension  coupling  matrix  and  the 

zero  k  and  k  curvature  terms.  When  any  of  B-,  •  or  Bfii  (i  =  1,2,6) 

terms  are  not  zero,  a  normal  strain  would  result  in  an  externally 

applied  M  or  M  moment  resultant  to  maintain  a  zero  curvature, 
x  xy 

In  this  paper,  we  will  be  concerned  with  symmetric  laminates, 
which  do  not  have  bending-extension  coupling.  If  we  consider  the 
specific  case  of  a  symmetric  laminate  under  uniform  axial  extension, 
ex  =  e0,  we  can  take  advantage  of  the  following  symmetry  conditions 
for  the  x-y  plane 

u (x,  y,  z)  =  u (x ,  y,  -z) 

v(x,  y,  z)  =  v (x ,  y,  -z)  (A-7) 

w(x,  y,  z)  =  w(x ,  y,  -z) 

and  for  the  x-z  plane 

v (x ,  y,  z)  =  -v(x,  -y,  z) 

(A- 8) 

w(x,  y,  z)  =  -w(x,  -y,  z). 

Because  of  the  above  symmetry  conditions,  a  symmetric  laminate  under 
uniform  axial  extension  requires  only  one  quarter  of  the  y-z  plane 
to  be  modeled. 

Since  a  linear  elastic  material  was  assumed,  the  constitutive  law 
can  be  expressed  as 


=  C^Uj  -  YT) 

i ,  j  =  1 ,  6 

(A-9) 

=  S^Oj  +  cj4T 

i,  J  =  1,  6 

(A- 10) 
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where  o-  =  stress  components 

e.  =  strain  components 
J 

C..=  modulus  matrix 
ij 

S..j=  compliance  matrix 

ct.j  =  coefficients  of  thermal  expansion 

AT  =  temperature  change 


From  equation  (A-lO)the  first  equation  (i  =  l)  can  be  rewritten  in  terms 
of  a-j  ,the  axial  normal  stress  component,  as 


al  ^e1  "  ATal  ‘  S12a2  “  S13a3  '  S16°6^ 


( A- 11 ) 


noting  is  the  axial  strain.  Since  a-j  can  be  expressed  in  terms 
of  the  axial  strain  and  the  other  five  stress  components,  then  equation 
(A-9)  can  be  expressed  in  a  reduced  form  as  [12], 


=  Dij(ej  -  ajAT)  +  at  i,  j  =  1,  5  (A-12) 

where  6^  =  o(i  +  1 ) 

Ej  =  e ( j  +  1 ) 

Ct j  =  a(  j  +  1 )  i  ,  j  =  1  >  5 

Did  =  C(i  +  l)(j  +  1) 
o*  =  (£0  '  aTAT)ci (j  +  1) 

With  the  reduced  constitutive  law  established  as  in  equation 
(A- 12),  the  finite  element  formulation  is  conventional.  Eight-noded, 
isoparametric,  quadrilateral  elements  were  used  having  conventional 
quadratic,  serendipity  shape  functions  [15]  and  three  degrees  of 
freedom  (u,  v,  and  w)  per  node.  The  global  system  of  simultaneous 
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linear  algebraic  equations  is 

K  ij-qj  =  QfT  +  Q-*  i,  j  =  1,  NDOF  (A-13) 

where  K.. .  =  global  stiffness  matrix 
q.  =  nodal  displacement  vector 

vJ 

=  in-plane  thermal  load  vector 
a* 

=  out-of-plane  load  vector 
NDOF  =  number  of  degrees  of  freedom 


The  above  system  of  equations  is  assembled  from  their  elemental  constit- 


uents  which  are 

K? .  =  f E 

U  vJ 

!  .D  B  .dV 
nn  mn  nj 

0e(AT)  _ 
wi  v 

/*B  .D  a  dV 
rj  mi  mn  n 

i,  j  =  1,  24 
m,  n  =  1 ,  5 

(A-14) 

Qe^*)  =  - 

/  B  -  D  a*  dV 
/  mi  mn  n 

where  B  .  =  the  nodal  displacement  strain  transformation 

•  I  j 

Dm  =  reduced  modulus  matrix  (equation  A-12) 

an  =  reduced  coefficient  of  thermal  expansion  (equation  A-12) 

★ 

a  =  in-plane  stresses  due  to  out-of-plane  strain  and  temperature 
change  (equation  A-12) 

V  =  the  volume  of  the  element 

The  boundary  conditions  for  the  problem  in  Figure  A3  can  be 
easily  established  by  symmetry;  these  are 
u(o,  z)  =  0 

v(o,  z)  =  0  (A-15) 

w(y,  o)  =  0 

This  allows  only  one  quarter  of  the  laminate  to  be  modeled. 


Figure  A3.  Finite  element  representati 
of  specimen  cross-section. 
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Verification  of  Finite  Element  Program 

Considerable  time  was  spent  verifying  the  finite  element  program 
due  to  conflicting  solutions  found  in  the  literature. 

Stress  distributions  near  a  traction  free  boundary  are  not  the 
prime  concerns  of  this  study,  but  the  solution  for  stresses  near  a 
traction  free  boundary  should  provide  a  good  check  for  the  program. 
Wang  and  Crossman  [12,30,31]  have  used  a  similar  uniform  axial  exten¬ 
sion  element  formulation  to  study  edge  effects  associated  with 
mechanical  and  thermal  loading  with  transient  moisture  sorption. 

The  first  problem  considered  was  the  stress  distributions  near  a 
free  edge  due  to  a  purely  mechanical  load  as  in  Reference  [12].  The 
following  two  stacking  sequences  were  considered;  [0/90]s  and  [90/0]s. 
The  laminate  geometry  is  noted  in  Figure  A4.  Due  to  the  symmetrical 
loading  and  stacking  sequences,  only  one  quarter  of  the  cross-section 
requires  modelling. 

The  load  is  input  as  a  unit  axial  strain  e  =  e  =  1 .  The  ther- 

X  0 

mal  load  is  set  to  zero.  The  ply  thickness,  h,  to  half  width,  b, 
ratio  was  1:8,  Figure  A4(b).  The  material  properties  used  were: 

Ex  =  137.8  GPa;  E2  =  E3  =  14.5  GPa 

G^  =  =  G22  =  5.86  GPa 

v12  =  v13  =  v23  =  0,21 

The  mesh  used  consisted  of  120  elements  with  405  nodes.  Each 
ply  contained  five  elements  through  the  thickness.  The  element 
aspect  ratio  was  kept  at  a  maximum  of  5:1. 

A  comparison  of  the  present  program  and  the  results  obtained  from 
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Reference  [12]  are  shown  for  a  stress  distribution  at  z  =  h  and  z  =  0 
in  Figures  A5  through  A10.  The  values  shown  from  Reference  [12] 
were  digitized  directly  off  of  a  photo  copy  of  the  article.  The  stress 
values  plotted  from  the  present  work  were  obtained  at  the  nearest 
Gaussian  integration  point  just  below  the  ply  interfaces  for  stresses 
at  z  =  h  and  just  above  the  midplane  for  stresses  at  z  =  0.  All  of 
the  figures  show  good  agreement  with  Reference  [12]. 

As  a  warning  to  the  reader,  caution  should  be  exercised  when 
examining  the  results  of  the  [±45]$  laminate  presented  in  Reference 
[12].  Some  of  the  plots  shown  were  drawn  with  an  incorrect  sign. 

The  readers  are  referred  to  Reference  [32]  where  the  correct  results 
for  the  free  edge  stresses  for  the  [±45]s  laminate  are  presented. 

The  next  loading  condition  was  a  thermal  load.  The  load  consis¬ 
ted  of  a  1°C  temperature  change  and  the  average  value  of  o 

XX 

transversely  across  the  laminate  was  zero.  The  constraint  of  the  zero 
average  value  of  cr  was  to  produce  a  purely  thermal  load  within  the 
confines  of  the  generalized  plane  strain  finite  element  program.  The 
same  material  properties  and  mesh  were  used  for  the  thermal  load  as 
were  used  for  the  mechanical  load.  Only  a  [0/90]$  stacking  sequence 
was  considered  for  the  thermal  load  case. 

A  comparison  of  the  present  work  and  the  results  obtained  by 
Reference  [30]  for  the  w  and  v  at  the  top  surface  and  free  edge, 
respectively,  are  presented  in  Figures  All  and  A12.  Stress  distribu¬ 
tions  at  the  midplane  (z=0)  are  shown  in  Figures  A13  through  A15.  The 
values  shown  for  Reference  [30]  were  obtained  in  the  same  manner 
as  stated  for  the  mechanical  loading.  Good  agreement  is  again  shown 
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Figure  A6.  a  distribution  at  z  =  h,  for  a  [90/0]  laminate 
under  mechanical  loading. 


W/h  10  u° 


Figure  All.  Displacement  W  at  z  =  2h,for  a  [0/90] 
laminate  under  thermal  loading. 


10 


.0 


(1  -  Y/b) 


Figure  A13.  a  distribution  at  z  =  0, in  a  [0/90] 

A  A 

laminate  under  thermal  loading. 


Figure  AI4*-cr  distribution  at  z  =  0,  in  [0/90] 
laminate  under  thermal  loading. 
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for  all  the  figures  in  the  thermal  load  case. 

Thus  far,  the  finite  element  program  has  demonstrated  its  ability 
to  predict  displacements  and  stresses  for  both  mechanical  and  thermal 
load  conditions.  The  remaining  area  to  be  considered  will  be  the 
fracture  analysis  capabilities  of  the  program.  This  area  will  be 
addressed  next. 

The  problem  considered  is  a  crack  in  an  isotropic  plate  of  finite 
height  and  infinite  width.  The  boundary  conditions  are  shown  in  Figure 
A16(a).  The  plane  solution  given  by  Rice  in  Reference  [33]  is: 


G 


II 


=  0 


This  problem  was  modeled  by  a  finite  width  plate  with  an  aspect  ratio 
of  b/h  =  5  as  shown  in  Figure  A16(b).  The  mesh  consisted  of  84  ele¬ 
ments  with  303  nodes.  The  crack  extension,  Aa  was  0.046b.  A  total 
of  15  crack  extensions  were  made.  The  results  are  shown  in  Figure  A17. 
The  value  of  Gj  is  shown  to  quickly  converge  to  Rice's  solution  at 
a  =  0.36b  and  remains  constant  for  the  remainder  of  the  solution. 

This  concludes  the  verification  of  the  finite  element  program. 
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Figure  A17.  Solution  to  a  crack  in  an  isotropic  plate 


of  finite  height. 


APPENDIX  B 


FATIGUE  TEST  DATA 


Fatigue  test  data  were  obtained  from  Ho  and  Schapery  [10].  The 
[±45/902^  AS/3501-6  samples  were  tested  at  24°C  and  50%  RH.  The  test 
frequency  was  3  Hz  with  a  load  ratio  of  +0.1.  Three  peak  fatigue 
stress  levels,  50%,  55%,  and  60%  of  the  ultimate  tensile  strength, 
were  used.  The  ultimate  tensile  stress  was  reported  as  135.6  MPa.  In 
Figures  Bl,  B2,  and  B3  the  mean  strain  curves  are  presented.  Several 
data  points  were  obtained  from  each  specimen  while  the  test  was  in 
progress.  Different  test  specimens  were  noted  by  the  different 
symbol s . 

The  solid  line  shown  is  the  mean  strain  function  used  in  the 
analytical  work  of  this  thesis.  The  function  was  presented  as  a 
piecewise  linear  curve  for  simplicity  and  to  avoid  any  bias  associated 
with  an  assumed  functional  form,  shown  together  in  Figure  B-4. 

Other  pertinent  data  relating  to  the  laminate  are  tabulated  in 
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Figure  B1 .  Mean  fatigue  strain  curve  for  a  [±45/902ls 
laminate  at  50%  UTS  maximum  fatigue  stress. 
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Figure  B2 .  Mean  fatigue  strain  curve  for  a  [±45/9C>2]s 
laminate  at  55%  UTS  maximum  fatigue  stress. 
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Figure  B3.  Mean  fatigue  strain  curve  for  a  [±45/902]s 
laminate  at  60%  UTS  maximum  fatigue  stress 
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Table  B1  AS/3501-6,  [±45/90J  .  Laminate  Data 
at  24°C  and  50%  RH  [10]  c  s 


Longitudinal  Young's 


Modulus,  E, 

«  1 

19.15 

GPa 

Ultimate  Tensile  Strength, 

UTS 

135.5 

MPa 

Ultimate  Strain 

10540 

m 
**  m 

Fiber  Volume 

(Process  Control  Panel) 

68.2% 

Fiber  Volume 
(Vendor  Data) 

62% 

Test  Specimen  Width 

20.63 

mm 

Ply  Thickness 

0.127 

mm 
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